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Background & Motivation 


Reacting Flows 


When chemical kinetic timescales are approximately equal to flow 
timescales, the chemical composition of a flowfield must be determined 
as part of a simulation procedure. Such flows are in chemical 
nonequilibrium. 


The physical models and governing equations for flows in 
thermochemical nonequilibrium have been simulated previously with 
finite difference and finite volume techniques. 

In this work we review the physical models and implement a SUPG finite 
element scheme for hypersonic flows in thermochemical nonequilibrium. 


T 
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capabilities of reusable thermal protection system materials. 

► Reusable materials typically limited to T < 2,000K. 

► It is necessary then to consider ablative materials. 

Ablative materials respond to high temperatures through pyrolysis, 
decomposition, blowing, and surface recession. 

Typically, ablation analysis is decoupled from the external flowfield, but 
we hope to do better. 

Additionally, accurately characterizing ground test facilities requires 
increased fidelity. 
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Governing Equations 


Governing Equations 

• Extension from a single-species calorically perfect gas to a reacting 
mixture of thermally perfect gases requires species conservation 
equations and additional energy transport mechanisms. 

f + v-(o»)=o 

^ + v • (puu) = —VP + V ■ r 
ot 

+ V • (pHu) = -V • q + V • (tm) 
ot 
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Governing Equations 
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V ■ ( pD s Vc s ) + lo s 
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^ + V • (puu) = -VP + V ■ r 
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dpE 


+ V • (pHu) = -V • q + V • (m) + V • ( p ^ h s V s Vc 
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Governing Equations 

• Extension from a single-species calorically perfect gas to a reacting 
mixture of thermally perfect gases requires species conservation 
equations and additional energy transport mechanisms. 

^ + V ■ (p s u) = V ■ ( pV s Vc s ) + u s 
at 

^ + V • (puu) = -VP + V ■ r 

Ot 

dpE ( ns 

+ V • {pHu) = -V ■ q + V • (tk) + V ■ ( p ^ h s V s Vc s 

\ i=i 


• Problem class may also require a multitemperature thermal 
nonequilibrium option. 


dpey 

dt 


+ V • ( pe v u ) = -V ■ q v + V ■ [ p ^ e Vs V s Vc s 


■ LOy 


5=1 
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Turbulence Modeling 

• We model the effects of turbulence using the Spalart-Allmaras 
one-equation turbulence model: 


^(/OZ'sa) + -^{pUjV sa ) =C b lS sa pU sa - C wi f w p 


d 


\d_ 

a dxk 


(p + pv m ) 


dv sa 

dx k 


Cb2 _ di2 sa dv$ a 

<7 dxk dxk 


with closure terms 

Pi = pVsnfvl , 


f w 

Jvl -j | 3 ’ 

X 3 + Si 


/v 2 = 1- 


X _ i'sa 

1 + Xfvl ’ X V ’ 


fw=g 


« 6 + c ii 


1/6 


g = r + c w2 r - r 


S sa K 2 d 2 
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Turbulence Modeling 

• We model the effects of turbulence using the Spalart-Allmaras 
one-equation turbulence model: 


9 {pv sa ) + sa ) =C hl S sll pV sll - C w \f w p 


dt 


I A 

a dxk 


/ . - \ Sr'sa 

[P + PP sa) 


Cb2 _ s a dVsi 


a ^ dx/c dxk 


and source term 
where 


S sa = n + s m , s m 0 = -^fvi 


S m = 


SmOj 

- ^ ^v3^mo) 

( (t’v3 2c v 2)tt 5/no) 


^ mO C v 2& 


otherwise. 
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Physical Modeling 


Thermochemistry 


Thermodynamics & Transport Properties 

• Thermochemistry models must be extended for a mixture of 
vibrationally and electronically excited thermally perfect gases. 

e int =e trans + e rot + e vib + e elec + h° 

ns 

=E^r s on+ E c seTm+ 

s— 1 s—mol 

ns ns 

E ^r b (T V ) + E c^r (Tv) + E c A° 

s—mol S— 1 5=1 

Here we have assumed that 7 n,ans = T ml = T and 7' vlb = 7’ elec = T v 

• Additional transport property models are required. In this work we use 

► species viscosity given by Blottner curve fits, 

► species conductivities determined from an Eucken relation, 

► mixture transport properties computed via Wilke’s mixing rule, and 

► mass diffusion currently treated by assuming constant Lewis number. 
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Physical Modeling 


Thermochemistry 


Chemical Kinetics & Energy Exchange 

Kinetics: 

• we consider r general reactions of the form 

N 2 + M — 2N + M 

N 2 + O ^ NO + N 

• When combined with forward and backward rates, these reactions 
produce the species source terms u s 

• Presently, we use either Canter A or an in-house library to provide 
these source terms. 

Energy Exchange: 

• Equilibration between the energy modes is modeled with a typical 
Landau-Teller vibrational energy exchange model with Millikan-White 
species relaxation times. 

• Provides the vibrational energy source term uy 

m 

September 26, 2012 11 /52 


Kirketal. (NASA/JSC) 


Fully Implicit Methods for Hypersonics 




Physical Modeling 


Thermochemistry 


Chemical Kinetics 

• We consider r general reactions of the form 

N 2 + M — 2N + M 


N 2 + O ^ NO + N 


• The reactions are of the form 


TZ r = k br ] 

.5=1 



kfr n 



a s 


where a sr and (3 sr are the stoichiometric coefficients for reactants and products 
• The source terms are then 


nr 

Us = Ms ^2 ( a ^ “ M (Kbr - 1 Zfr) 
r=l 
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Physical Modeling 


Thermochemistry 


Energy Exchange 


& V — Qv H - (^transfer 
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Physical Modeling 


Thermochemistry 

Energy Exchange 

tU V — Qv T ^transfer 

We adopt the Landau-Teller vibrational energy exchange model 


er vib = p,— 

1 s 

where ej lb is the species equilibrium vibrational energy and the vibrational 
relaxation time r s vlb is given by Millikan and White 


T. Vib = 


E ns 

r=l Xr 

Y^UXrjrl 


M I c s 

= c ' mT’ m = £ w. 


and 


^ exp [A sr (t" 1 / 3 - 0.015m!/ 4 ) - 18.42] 


A sr = 1.16 x 10 V!/X/ 3 , Msr = 


M s M r 
M s -f- M r 
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Physical Modeling 


Thermochemistry 


Vibrational Energy Production and Energy Exchange 

— Qv + Q transfer 
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Physical Modeling 


Thermochemistry 


Vibrational Energy Production and Energy Exchange 

&V = Q v + <2 transfer 

When molecular species are created in the gas at rate 6j s , they contribute 
vibrational/electronic energy at the rate 

Q v ,=*,(e?> + e?*) 

so the net vibrational energy production rate is 

ns 

e,. = X>K ib + ^ lec ) 

5=1 
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Physical Modeling 


Thermochemistry 


Vibrational Energy Production and Energy Exchange 

&V = Q v + <2 transfer 

When molecular species are created in the gas at rate 6 j s , they contribute 
vibrational/electronic energy at the rate 

G»=^(^ b + ef c ) 

so the net vibrational energy production rate is 

ns 

e,. = X>K ib + ^ lec ) 

Combining terms yields the desired net vibrational energy source term 

ns ns 

^v = Eer ib + E^ vib +ef c ) 

5=1 5=1 
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Physical Modeling 


Quasi-Steady Ablation 


Ablation Processes 


Pyrolysis Gas 



Ablating 
/ Surface 


iiL 

A-' • v.v ‘7' ..V .. 'A 
Sv.-. :• vV 



. Virgin Composite » * 
/Sub-Structure 

J_LLLL 

Schematic of ablation processes 


Char Zone 


- Decomposition 
Zone 


• Ablation is a multi-scale, multi-physics phenomenon 

• Sometimes amenable to simplification for predictive simulations 
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Physical Modeling 


Quasi-Steady Ablation 

Quasi-steady State Ablation Hypothesis 



O Steady state in reference frame fixed to the receding surface 
© Time variations solely due to motion of the material domain 

0 Time scale for surface motion (s ~ 0. 1 — I mm / sec ) much larger than 
characteristic time scale of unsteady processes 

Q 1-D, semi-infinite medium 
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Quasi- steady Ablation 


Quasi-Steady Ablation 


Energy 

Mass 



p,y a 


Pv V cshf,v 


Zone Material 


• Assumes ablation timescale <C trajectory timescale 

• Assumes negligible substructure conduction 
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Physical Modeling 


Quasi-Steady Ablation 

Ablation Interface Conditions 

Recession: 

pv w = m c + m g 

Mass: 

Ji\ gas + pv w Cj = Nj(Cj, T ) + m g C itg ; ( i : 1 ..N s ) 

Energy: 

dT I Ns „ 

~ k-y: ~ ^2 hj(T w ) Jj\ gas + m'h c (T ) - pv w h w {T ) 

r l g ai j = i 

% „ 

-\-OiCj r — <JtT w 4 + ^ ^ m g Ci >g hi(T w ) + 

o I solid., xv — 0 
i=l ^ 

• Nonlinear Robin Boundary Conditions 

• Enables quasi-steady solves, restarts 
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Stabilized Finite Element Scheme 


dU dFj dGj ■ 

1 = — - + S 

dt dxj dxj 
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Stabilized Finite Element Scheme 


dU dU 

dt 1 dxi 


d_ 

dxi 



+ S 
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Finite Element Formulation 


Stabilized Finite Element Scheme 


dU dU 

dt 1 dxi 


d_ 

dxi 



+ S 


Find U satisfying the essential boundary and initial conditions such that 



+ 


dW 

d X: 


dU 

K ^r A ' v 


'dU dU dGi 

dt 1 dxj dxi 


dn 

dtt 



l 


W ■ ( g - f ) 


dT = 0 


for all W in an appropriate function space. 
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Finite Element Formulation 


Stabilization Parameters 

Discontinuity capturing operator: 


D3CO = 


au , a du a ( K au\ 

at t; dxi y^iidxj) 

2 

V 


4 -l 0U h 
^0 dxj 

SUPG stabilization matrix: 



t supg — 


£ 

*=nodes 


d<l>i A d<t>j \ 

dxj 7 dxj ik dxu ) 


where 



= L\A\R 
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Finite Element Formulation 


Fully-Implicit Navier-Stokes (FIN-S) 


Implementation Highlights 

• C++ application code built on top of the libMesh library. 

► libMesh provides all requisite finite element data, parallel domain 
decomposition details. 

► Inherits PETSc preconditioned Krylov iterative solvers. 

► Cantera used for kinetic rates, in-house thermodynamics, transport 
properties. 

► Only ss 30K SLOC 

• Fully-coupled (monolithic solves), fully-implicit discretization. 

• Rigorous verification using MASA-provided manufactured solutions. 

• Testbed for intrusive VV/UQ schemes applied to hypersonics. 


September 26, 2012 22/52 


Kirketal. (NASA/JSC) 


Fully Implicit Methods for Hypersonics 




Finite Element Formulation 


Spalart-Allmaras Perfect-Gas Verification 
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Parallelism 


Need for Parallelism 


Large Problem Size 

• Large numbers of unknowns. 

► For a Lagrange nodal basis: 

# DOFS = (NS + NDIM + NE + NT) x # NODES 

► Specifically, for our 13 species ablation model in 2D with turbulence 

# DOFS = (13 + 2 + 2+ 1) x# NODES 

• For our implicit scheme, both storage and computational cost scale like 
(# DOFS) 2 
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Parallelism 


Need for Parallelism 


Complex Physical Models 

• Chemical Kinetics, transport properties for NS species inherently 
expensive. 

• Temperature is a nonlinear function of species concentration, internal 
energy for a mixture of thermally perfect gases. 

• Quasi-steady ablation boundary condition is also nontrivial. 
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Parallelism 


Opportunities for Parallelism 

Multiple Types of Parallelism 

O Domain Decomposition: We use a standard non-overlapping domain 

decomposition approach provided by libMesh. Local computations are perfectly 
parallel, and the resulting implicit system is solved using preconditioned Krylov solvers 
from PETSc. 

© Multithreaded Computation: The relatively large element matrices resulting for 
reacting flows are well suited for threaded assembly. libMesh provides a convenient 
interface to Intel’s Threading Building Blocks which can provide further parallelization 
on multicore architectures. 

© Vectorization: Remember vectorization? While no longer the de facto paradigm for 
high-performance computing, modern microprocessors offer vectorized instructions 
worth exploiting. We are using Eigen for dense linear algebra and inherit its SSE 
optimizations. 
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Domain Decomposition 


Parallelism 
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Domain Decomposition 



ids for Hypersonics 
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Parallelism 


Speedup - Domain Decomposition 





Multithreading 


Parallelism 


• Modem Parallel systems often contain 12-16 (or more) on-node cores 
connected via low-latency network. 


• On-node multithreading allows an additional parallel mechanism that 
can extend scalability in certain circumstances. 

• libMesh provides a clean interface to Intel®’s Threading Building 
Blocks (TBB) which is we have access to. 


• TBB is a C++ template library consisting of 

► Algorithms 

► Containers 

► Mutexes 

► Timing routines 

► Memory allocators 


designed to help avoid low-level use of platform-specific (e.g. 
pthread) implementations. 
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Parallelism 


Intel® ’s Threading Building Blocks 


• Requires more work than OpenMP but 

► Has better type-safety 

► Easier to reuse code 

► More natural for use with C++ 

• Once a standard for loop is selected for parallelization its components 
are abstracted as C++ Range and Body objects 


In FIN-S we parallelize matrix assembly, primitive variable computation, 
and other operations in this way. 

► Some operations perfectly asynchronous - e.g. computing primitive 
variables. 

► Other operations require locking shared objects - e.g. inserting local 
contributions to a global matrix. 

► Special care needed when interfacing with 3 rd party libraries. 
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Results 


Viscous Thermal Equilibrium Chemical Reacting Row 

2D Extended Cylinder 

• Laminar flow in thermal equilibrium 

• Chemical nonequilibrium, 5 species air (N 2 , CL, NO, N, O) 

• 5 reaction model with Park 1990 rates 

cN2 )00 = 0.78,cO2,oo = 0-22 

Coo = 6 , 7 3 1 m /sec 
Poo = 6.81 X 10 ^ kg/m 3 
Too = 265 K 

• Blottner/Wilke/Eucken with constant Lewis number Le = 1.4 for 
transport properties 

• Mesh, iterative convergence 

• FIN-S/DPLR comparison 
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Results 


Viscous Thermal Equilibrium Chemical Reacting Row 




Nitrogen 
Mass Fraction 
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y (rn) 


Results 


Viscous Thermal Equilibrium Chemical Reacting Row 


Code-to-Code Comparison - 

Stagnation Line 




x(m) 



Stagnation Line 
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Results 


Viscous Thermal Equilibrium Chemical Reacting Row 


Code-to-Code Comparison - 

Flank Line 
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P (N/m 2 ) 


Results 


Viscous Thermal Equilibrium Chemical Reacting Row 


Mesh Convergence 






400x400 

200x200 

100x100 


Kirketal. (NASA/JSC) 


Fully Implicit Methods for Hypersonics 


x(m) 

September 26, 2012 


40/52 




Relative Transient Residual, IAC//A/I 


Results 


Viscous Thermal Equilibrium Chemical Reacting Row 


Iterative Convergence 
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Time Step Size, At, (seconds) 



Results 


Viscous Reacting Flow with Quasi-Steady Surface Ablation 

Ablating Boundary Experiments 



• Turbulent flow in thermochemical nonequilibrium, 13 species air (N 2 , 

0 2 , NO, N, O, C 3 , C 2 , C, CN, CO, H 2 , H, C 2 H), 18 reaction model with 
Park 2001 rates 

• 5 Meter-scale domain, millimeter-scale chemical boundary layer 
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Results 


Viscous Reacting Flow with Quasi-Steady Surface Ablation 


Ablating Boundary Experiments 





Arcjet Flowfields 


Results 


Modeling Arcjet Flows 


Motivation 

• Arcjets are uniquely suited to perform high enthalpy, long duration 
material response testing. 

• Modern computational techniques are required to adequately 
characterize the freestream properties. 

• Analysis complicated by multitude of scales, physical phenomenon: 

► Very low speed, high pressure plenum, 

► very high speed, low pressure nozzle exit, 

► highly nonequilibrium flow about test specimen. 

• Adequately treating these phenomenon simultaneously is challenging for 
numerical methods. 
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Arcjet Flowfields 


Results 


Modeling Arcjet Flows 
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Arcjet Flowfields 


Modeling Arcjet Flows 






Arcjet Flowfields 


Results 


T(K) 


JSC TP2 
5" nozzle 
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Arcjet Flowfields 


H 0 (J/kg) 


JSC TP2 
5" nozzle 


Modeling Arcjet Flows 




s for Hypersonics 
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Ongoing Challenges 


Full Disclosure 


Opportunities for Further Enhancement 

O Linear Solver Strategy: Preconditioned GMRES is highly effective but 
potentially overkill for early, highly nonlinear transients. Mixed implicit/explicit 
schemes may provide a fast alternative. 

0 Improved Shock Capturing: Robust shock capturing is still a challenge. Current 
scheme is fragile on bad meshes, and often convergence stalls. 
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Additional Focus Areas 
O Physics Modeling 

► Weakly Ionized Flows 

► Additional turbulence models 

► Fully coupled radiative transport 

0 Unsteady ablation coupling 
0 Adjoints 

► Sensitivity analysis 

► Adaptivity 




Ongoing Challenges 


Thank you ! 
Questions? 
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